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ABSTRACT 

Aims. We aim to demonstrate the efficiency of a Bayesian approach in analysing radial velocity data by reanalysing a 
set of radial velocity measurements. 

Methods. We present Bayesian analysis of a recently published set of radial velocity measurements known to contain 
the signal of one extrasolar planetary candidate, namely, HD 11506. The analysis is conducted using the Markov chain 
Monte Carlo method and the resulting distributions of orbital parameters are tested by performing direct integration 
of randomly selected samples with the Bulirsch-Stoer method. The magnitude of the stellar radial velocity variability, 
known as jitter, is treated as a free parameter with no assumptions about its magnitude. 

Results. We show that the orbital parameters of the planet known to be present in the data correspond to a different 
solution when the jitter is allowed to be a free parameter. We also show evidence of an additional candidate, a 0.8 Mjup 
planet with period of about 0.5 yr in orbit around HD 11506. This second planet is inferred to be present with a high 
level of confidence. 
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1. Introduction 

The question of whether an extrasolar planet is detectable 
or not depends on a delicate mixture of observational tech- 
nology, effective tools for data analysis, and theoretical un- 
derstanding on the related phenomena. Traditionally, the 
instrumentation has been the most celebrated part of this 



trinity (Santos et al. 2004 Li et al. 20081, but the fainter 



the signals that one is able to detect, the more the attention 
should be directed towards the two other factors involved 
in successful discoveries. It is already known that of poor 
understanding of the noise-generating physics may produce 
misleading results ( Queloz et aH|20C)l| , although this is also 
true for the means the data analysis. In many cases the sta- 
tistical methods involved in the data reduction and analysis 
of the observations have a tendency to identify local solu- 
tions instead of global ones. This is particularly true espe- 
cially when information is extracted using gradient-based 
algorithms. 

In this letter, we reanalyze the radial velocity (RV) data 
of a detected extrasolar planet host, namely HD 11506 
(Fischer et al. 2007| ). Our analysis is based on Bayesian 



model probabilities and full inverse solutions, i.e., the full 
a posteriori probability densities of model parameters. The 
Bayesian model probabilities provide a strict mathemati- 
cal criterion for deciding how many planetary signals are 
present in the data. We calculate the inverse solutions us- 
ing the Markov chain Monte Carlo (MCMC) method with 
the Metropolis-Hastings algorithm (Metropolis et al. 1953 



iHastingsl 119701). These inverse solutions are used to present 



reliable error estimates for all the model parameters in the 
form of Bayesian confidence sets. As a result, we present 
strong evidence for a new and previously unknown plane- 
tary companion orbiting HD 11506. We also demonstrate 
the importance of treating the magnitude of stellar RV 
noise, the RV jitter, as a free parameter in the model. 



2. The model and Bayesian model comparison 

When assuming the gravitational interactions between the 
planetary companions to be negligible, the RV measure- 
ments with k such companions can be modelled as (e.g. 
GreenI flgSSl) 



m 



COs(l^j(t) -|- UJ-j) + Cj COS LOj 



(1) 
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where Vj is the true anomaly, Kj is the RV semi-amplitude, 
Ej is the eccentricity, and ujj is the longitude of pericentre. 
Index j refers to the jth planetary companion. Hence, the 
RV signal of the jth companion is fully described using five 
parameters, Kj, Wj, e_,-, mean anomaly Mqj, and the orbital 

period Pj. 

Following Gregory (2005a 2007a|b l, we calculate the 
Bayesian model probabilities for different statistical mod- 
els. These models include the first and second companion 
Keplerian signals and an additive RV jitter (models Ml 
and M2, respectively). This jitter is assumed to be Gaussian 
noise with a zero mean and its deviation, aj, is another free 
parameter in the model. The jitter cannot be expected to 
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provide an accurate description of the RV noise caused by 
the stellar surface. It also includes the possible signatures of 
additional planetary companions if their signals cannot be 
extracted from the measurements. Hence, this parameter 
represents the upper limit to the true RV variations at the 
stellar surface. We also test a model without companions 
(MO) for comparison. 

The instrument error, whose magnitude is usually 
known reasonably well, is also included as a Gaussian ran- 
dom variable in our analysis. Hence, for a model with k 
planetary companions, the RV measurements at U are 
described by 



ri = Zk{U) + j + ei + ej, 



(2) 



where e/ ~ N{0, aj) is the instrument error, ej ^ N{0, ctj) 
describes the remaining uncertainty in the measurements, 
called the RV jitter, and parameter 7 is a reference velocity 
parameter. 

When deciding whether a signal of a planetary compan- 
ion has been detec ted or not, we adopt the Bayesian model 
selection criterion ( Jeffreys 1961 1 . This criterion states that 
a model with k + 1 planetary companions should be used 
instead of one with k companions if 



P(ifc+i|r) » P(ife |r), 



(3) 



where ris a vector consisting of the RV measurements. For 
more informatio n ab out Bayesian m odel compa rison, see 



Kass fc Raftery (e.g. |1995[ ); [Gregory] (e.g. |2005a l. 

After identifying the most appropriate model in the con- 
text of Eq. ([3]) , we must find the Bayesian credibility sets 
that can assess accurately the uncertainties in the orbital 
parameters and the RV masses of the planetary compan- 
ions. The Bayesian credibility sets are robust uncertainty 
estimates because they show the uncertainty of the model 
parameters given the measurements. A Bayesian credibil- 
ity set Vs containing a fraction S G [0, 1] of the probability 
density of parameter vector u can be defined to a subset 
of the parameter space U that satisfies the criteria (e.g. 
io & Somersalol 



[Kaipic 



20051 



uevs Piu\m)du = 
p{u\m)\ueavs = c. 



(4) 



where p(u\m) is the probability density of the parameters, 
c is a constant, m is the measurements, and dVg is the edge 
of the set Vs. We use S = 0.99 throughout this article when 
discussing the parameter errors. 

The question of model comparison is not only statistical 
but also physical. Because the inverse solution is based on 
non-interacting planets - a simplification that is fairly ad- 
equate in most cases - it might be physically impossible in 
reality. This is important to understand because many ex- 
trasolar systems contain large eccentricities that may cause 
close encounters. In these detailed description of 

the dynamics is required. We drew 50 random values of u 
from its posterior probability density and integrated these 
directly to investigate the orbital evolution. The inte gration 
was conducted using the Bulirsch-Stoer integrator (Bulirsch 



& Stoer 1966), which is famous for its high reliability even 
when the dynamics include several close encounters. The 
stability analysis that we used in this analysis is not exact 
but estimates reliably the expected long-term behaviour 
from a relatively short ensemble or orbit. In this simple 
analysis, we identify the variation in the semi-major axis 



and the eccentricity from randomly selected initial condi- 
tions to the end of integration. If the variations are not sig- 
nificant the system is then assumed to be long-term stable 
(N. Haghighipour, priv. comm.). The advantage this formu- 
lation is the short integration time required when compared 
to more adequate, for example, Lyapunov exponent based, 
integration methods. 

The most important error source in this type of study 
is a phenomenon known as stellar jitter. It is caused by a 
combination of convection, rotation, and magnetic activity 
on the stellar surface (see e.g., Wright, 2005, and references 
therein). Although the role of stellar activity as an error 
source is well known, the magnitude of this error is almost 
always assume d to be constant, based on certain studies 
(Wright 20051. Here, we use a more conservative approach 
and consider the magnitude of the jitter to be free param- 
eter. 

Before this work, four exoplanets were found using 
Bayesian approach instead of a more traditio nal peri- 
odogram, the first of them being HD 73526 c (Gregory 
2005 a[). The ca ndidate HD 208487 c was proposed by 



Gregory (2005b I and later confirmed by ad ditional observa- 



tions (Butler et al. 2006 



2007a I. For the system 



Gregory 

HD 11964, the Bayesian analysis revealed evidenc e of three 
planets, despite only one being previously known (I Gregory] 
2007b| ) . None of these works include a discussion of the dy^ 
namics, although it is clear that including dynamics in the 
work enhance the quality of results. This point also is inter- 
esting because if the dynamical analysis excludes a part of 
the parameter space as physically impossible, this restric- 
tion can be inserted into Bayesian model as additional a 
priori constraint that will, with the data, provide tighter 
confidence limits. 



3. Orbital solutions 

The star HD 11506 is a quiescent main sequence star of 
spectral class GO V. It is a relatively nearby star with a 
Hipparcos parallax 18.58 mas, which corresponds to a dis- 
tance of 53.8 pc. It has T ^ff = 6060 K and [Fe/H] = 0.31 
dValenti fc Fischer] 



^ , ,2005| and its mass is estimated to be 

i.l9 Me ( [Fischer et al.H2007t . We use this mass estimate 



throughout the paper when calculating planetary masses 

The pla netary companion H D 11506 b was first 
nounced by Fischer et al. ( 2007 1 . They speculated an addi- 



tional companion could be present because the value of 
their single-companion model fit was large (10.3). However, 
their value was c alcula ted by assuming a fixed jitter 



level. Fischer et al. (20071 assumed that aj = 2.0ms ^. 
Our solution for this parameter is consistent with this esti- 
mate (Table [2|). 

The Bayesian model probabilities for k — 1,2 are listed 
in Table [T] These probabilities imply that M2 provides tha 
most accurate description of the RV data of HD 11506. 
These probabilities favour two planetary companions, im- 
plying that, according to the measurements of [Fischer et[ 
al. (2007), this star does host two companions. It also ver- 
ifies that the large value was consistent with there be- 
ing a signal o f an ad ditional companion, as suspected by 



Fischer et al. 



([2007[). The corresponding x value of our 
two-companion solution is 3.1. They also mentioned that a 
170 day period could be present in the data but were unable 
to verify the existence of a second companion. This period 
corresponds to our two-companion solution (Table [2 ) . 
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Table 1. Bayesian model probabilities of 1 and 2 planet 
models. 



Model 


Probability 




HD 11506 


Ml 


< 10"" 


M2 


1 



Table 2. The RV two-planet solution of HD 11506. MAP 
estimates of the parameters and their I?o.99 sets. 

Parameter MAP O0.99 

Pi [yr] OS [3.22, 4.01] 

ei 0.22 [0.10, 0.47] 

Ki [rns"^] 57.4 [49.7, 71.2] 

uji [rad] 4.5 [3.8, 5.1] 

Ml [rad] 4.7 [2.5, 0.5] 

mp,isinji [MjJ 3.44 [2.97,4.34] 

ai [AU] 2.43 [2.31, 2.67] 



P2 [yr] 0.467 

62 0.42 

K2 [ms-i] 25.5 
LU2 [rad] 4.1 
Ma [rad] 5.5 

mp,2 sin 42 [Mj„p] 0.82 

a2 [AU] 0.639 



[0.450, 0.476] 
[0, 0.62] 
[11.2, 35.3] 
[2.3, 5.1] 
[0, 27r] 
[0.32, 1.13] 
[0.622, 0.646] 



7 [ms \ 
aj [ms"^] 



6 

3.5 



[-4, 15] 
[0.6, 8.8] 



W hen as suming the fixed jitter level adopted by Fischer| 



et al. (20071, the probability of the one-companion model 
was again found to be considerably lower than in Table [T] 
(less than 10"'^^). This result emphasizes the fact that the 
jitter cannot be fixed but must be considered as a free pa- 
rameter. Furthermore, by fixing the jitter to some a priori 
estimated value, the probability densities of the orbital pa- 
rameters become far narrower, which underestimates the 
uncertainty in the solution. 

On Fig. [1] we show the maximum a posteriori (MAP) 
orbital solution of model M2 and the velocity curve of HD 
11506 c after the signal of b companion has been subtracted. 
Interestingly, our s olution for the companion HD 11506 b 
differs from that of Fischer et al.| (^007). For instance, we 
found its RV mass to be 'SAWljup, whereas they reported 
a mass of 4.7Mj„p, which is outside the margins of I'd. 99 
set in Table |2] The jitter parameter also appears to have 
a higher value than estimated. This could be indicative of 
the difficulties in estimating the jitter magnitude but the 
jitter may also contain the signal of a third companion. 

We were unable to determine any strong correlations 
between the model parameters, and all of the parameter 
densities, apart from P2 (Fig. [2]), were reasonably close to 
Gaussian. 

We subtracted the MAP one-planet solution from the 
RV meas urements and calculated the Scargle-Lomb peri- 
odogram (Scargle 1982 Lomb 19761 for these residuals. 



The highest peak corresponded to a period of 0.45 yr, which 
is close to the period of the second companion (0.47 yr). 
However, the FAP of this peak was as high as 0.54, which 
ensures that it was impossible to detect the signal of this 
companion with periodogram. In contrast, this solution was 
easily found using the MCMC method. In agreement with 
the MAP solution, we were unable to find other probabil- 
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Fig. 1. Radial velocity measurements of HD11506 (Fischer 
et al. 2008). Two companion MAP orbital solution (top), 
and signal of companion HD 11506 c (bottom). 



ity maxima in the parameter space of the two-companion 
model. 



4. Orbital stability 

To present additional evidence of HD 11506 c, we selected 
50 random combinations of parameters taken from the pa- 
rameter probability densities. These values were used to 
simulate the dynamical behaviour of the planetary system 
by direct integration. A random example of these simula- 
tions is shown in Fig. |3] (top), where 100 000 yr excerpts 
are presented for two randomly selected parameter com- 
binations. The orbital ellipses precess slowly but both the 
semimajor axis and the eccentricities remain almost con- 
stant during the evolution. This feature does not prove but 
suggests strongly that the two planetary companions orbit- 
ing HD 11506 provide a physically stable system. 

We tested the stability further by selecting 50 values of 
the parameter vector from the region of parameter space 
most likely prone to instability, i.e., where the quantity 
Q = ai(l — 61)7(02(1 -1-62)) has its smallest value. We drew 
the values from a probability density, whose maximum is 
at minQ and that decreases linearly to zero at maxQ. This 
results in a sample from the posterior density that has val- 
ues mainly close to minQ. This test was designed to ob- 
tain further information about the orbital parameters by 
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Radidal Velocity of HD1 1506 




0.44 



0.46 



0.48 



Fig. 2. The posterior probability density of parameter P2 
and its mode, mean, deviation, and two higher moments. 
The sohd curve is a Gaussian curve with the mean and 
variance of the density. 




Fig. 3. Two random examples of two-planet solution for 
HD 11506. Green and blue areas indicate orbits of HD 
11506 b and c, respectively. Red cross denotes the star. 
Both figures present a 100 000 year part of orbital evolu- 
tion, and the slow precession of the apsid line is clearly 
visible. Unit of both axes is AU. 



excluding parameter values that excluded an unstable sys- 
tem. However, all 50 parameter values resulted in bound 
orbits even though they precessed strongly. Therefore, it 
was impossible to extract additional information about the 
orbital parameters from these simulations. See random ex- 
ample in Fig. [3] (bottom). 

5. Conclusions 

We have presented complete Bayesian reanalysis of radial 
velocities of HD 11506, and identified the orbital parame- 
ters of a previously unknown exoplanet candidate. 



These analyses demonstrated the importance of taking 
stellar jitter, the most important error source in RV mea- 
surements, into account as a free parameter in the analysis. 
As an unknown noise parameter, jitter cannot be predeter- 
mined to some estimated value because of its strong effect 
on the Bayesian model probabilities of models with different 
numbers of companions. If overestimated, it may prevent 
the detection of a additional companions. Its uncertainty 
must also be taken into account when calculating the error 
bars for the orbital parameters to prevent the underestima- 
tion of their errors. 

In the case of HD 11506 the two-pla net model is by 



far more probable for the given data (Fischer et al. 2007 1 



than the original one-planet solution, a result that remained 
unchanged regardless of whether we considered the jitter as 
fixed or free. The fit of the two-planet model is similar to 
any known two-planet signal and a dynamical analysis has 
demonstrated the solution to be physically possible. We 
therefore claim that the RV data, which was presented by 
|Fischer et al. (20071, contains the signals of two planetary 



companions. 

Another interesting feature of the full inverse solution of 
the two-companion model is that the RV mass of HD 11506 
b diffe rs significantly from that obtained by [Fischer et al. 
(2007). This difference implies that the Bayesian model se- 



lection criterion should be used when assessing the number 
of planetary signals in RV data. Without accurate knowl- 
edge on the best-fit model, the orbital solution can be 
biased and the resulting statistical conclusions about the 
properties of extrasolar planetary systems can be mislead- 
ing. 

More observations are required to tighten further the 
constraints on the parameter space of HD 11506 system. 
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